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ABSTRACT 

The first year Wilkinson Microwave Anisotropy Probe (WMAP) set quantitative con- 
straints on the amplitude of any primordial non-Gaussianity. We run a series of dark 
matter-only N-body simulations with the WMAP constraints to investigate the effect 
of the presence of primordial non-Gaussianity on large scale structures. The model 
parameters can be constrained using the observations of protoclusters associated with 
Ly-a emitters at high redshift (2 ^ z ^ 4), assuming the galaxy velocity bias can be 
modelled properly. High redshift structure formation potentially provides a more pow- 
erful test of possible primordial non-Gaussianity than does the CMB, albeit on smaller 
scales. Another constraint is given by the local galaxy density probability distribution 
function (PDF), as mapped by the 2 degree Field Galaxy Redshift Survey (2dFGRS). 
The PDF of 2dFGRS L^_^ galaxies is substantially higher than the standard model pre- 
dictions and requires either a non-negligible bias between galaxy and dark matter on 
~ 12 h~^ Mpc scales or a stronger non-Gaussianity than allowed by the WMAP year 
one data. The latter interpretation is preferred since second-order bias corrections are 
negative. With a lower normalisation of the power spectrum fluctuations, as = 0.74, 
as favoured by the WMAP 3 year data, the discrepancy between the Gaussian model 
and the data is even larger. 
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1 INTRODUCTION 

The best limits on possible initial non-Gaussianities in the 
primordial density fluctuation distribution have hitherto 
been imposed by the CMB (Komatsu & Spergel 2001). 
There are many models which can generate such non- 
Gaussianities, including multiple field inflation or more spec- 
ulative models such as those assuming non-homogeneous re- 
heating processes or those seeded with topological defects. 
Moreover such non-Gaussianities may be scale-dependent, 
and it is therefore important to complement CMB searches 
with smaller scale probes, on cluster and galaxy scales. 

Our approach here is phenomenological, following that 
of Mathis, Diego & Silk (2004), and focuses on the large-scale 
galaxy distribution. Tracers of primordial non-Gaussianities 
include possible large-scale power excesses in the nearby 
galaxy distribution as mapped by 2dFGRS and SDSS red- 
shift surveys, the early formation of protoclusters up to 
a ~ 4 as traced by radio galaxies surrounded by Ly-a emit- 
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ters and Lyman-break galaxies, and the possible excess of 
angular power in CMB temperature fluctuations seen by 
CBI and ACBAR on cluster scales (which is most likely 
due to the thermal SZ efi'ect). With the claimed reduc- 
tion in normalisation of the fluctuation power spectrum to 
(T8 = 0.74 ± 0.05 (Spergel et al. 2006), such questions take 
on renewed urgency. 



Our approach here is to study the implications via dark 
matter N-body simulations of a simple non-Gaussian model, 
which modifies the frequency of rare peaks in the primor- 
dial density field that defines the locations where massive 
clusters (or superstructures) are formed. If the probability 
distribution function of the density field shows an excess (or 
deficit) of high peaks compared to the Gaussian case, pro- 
toclusters of given mass form earlier (or later), and systems 
above a given mass will be more (or less) abundant at fixed 
redshift. Sadeh, Rephaeli & Silk (2006) presented analyti- 
cal models, whereas here we simulate and evolve a suite of 
N-body simulations of 300 h~^ Mpc on the side. 
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2 THE CASE FOR NON-GAUSSIANITY 

We consider a simple non-Gaussian model, where the am- 
plitude of non-Gaussianity in the primordial fluctuations is 
parametrised by a non-linear coupling parameter, /jv_t (Ko- 
matsu et al. 2003). We refer to this non-Gaussian model as 
the /jvi model. The Bardeen curvature perturbations $ at 
any given position x can be written as 



^X) = <i>L{x) + fNL[<i>Ux)~ < $1{X) >] 



(1) 



where $_L(a;) are Gaussian linear perturbations with zero 
mean. This parameterisation is useful for setting quantita- 
tive constraints on the amount of non-Gaussianity allowed 
by the CMB data. Komatsu et al. (2003) used first year 
WMAP data (WMAP-1) to establish limits on /nl and ob- 
tained —58 < /jvi < 134 at 95% confidence. 

We produce the initial conditions for the /nl simu- 
lations using the conventional grid method, but when cal- 
culating the displacement using the Zel'dovich approxima- 
tion, we use the perturbed newtonian potential, which is —1 
times the curvature perturbation given by Equation [1] We 
then run a series of dark matter N-body simulations using 
the publicly available code GADGET (Springel, Yoshida & 
White 2001). We consider three values for /nl (-58,0,134) 
consistent with WMAPl constraints, where /jvl = cor- 
responds to the standard Gaussian model. For these simu- 
lations, we impose the initial power spectrum to be given 
by an adiabatic CDM transfer function with spectral slope 
Us — 1. We use the 'old' concordance cosmological model. 
I.e. Qm = 0.3, r^A = 0.7, h = 0.7 and as = 0.9. Re- 
cent results from WMAP-3 prefer a lower value for erg 
(Spergel et al. 2006). We therefore run another series of sim- 
ulations using the WMAP-3 cosmological parameters (ie. 
fi™ = 0.238, = 0.762, h = 0.73, as = 0.74), both 
for /jvl = and /jvl ~ 134. All N-body simulations are 
started at redshift Zinu = 70 and we follow 128^ particles 
in a cube of 300 h~^ Mpc aside, implying a particle mass 
of ~ lO^^h~^M0, which, for the present study, gives high 
enough resolution as shown in Section |3] For each of the 
above models, we present the mean results from ten statisti- 
cally equivalent realizations. Hereinafter, we conservatively 
discuss the results for the cosmological model with erg = 0.9 
unless otherwise stated. 



3 RESULTS 

3.1 Halo mass function in non-Gaussian models 

As with any simulation run, it is essential to check if our 
analysis is affected by the chosen mass resolution. In Fig. [T] 
we show a comparison between our simulations and results 
from either other high-resolution N-body simulations or an- 
alytical models. The upper panel compares the halo mass 
function of our Gaussian simulations with the fitting for- 
mula of Jenkins et al. (2001), which is based on a series of 
high-resolution N-body simulations with similar cosmology. 
Our estimate of the halo mass function is in good enough 
agreement with the Jenkins et al. (2001) fit. In the lower 
panel, we show the probability distribution function (PDF) 
of the density contrast smoothed with a spherical top-hat 
filter of radius 12 h~^ Mpc. Kayo et al. (2001) have shown 



c 
-o 



-5 - 



-7 - 



This Simulation 
Jenkins el al. (2001) 



10'^ 

M (h-'Mg) 



10 w 



0.01 




1+6 

Figure 1. Upper panel: mass function of dark matter halos 
at z=0 from simulations with Gaussian initial conditions. The 
solid line shows the mean mass function from 10 simulations 
with Gaussian initial conditions, and the dashed line shows 
the Jenkins et al. (2001) halo mass function. Lower panel: 
the probability distribution function of the density contrast. 
Again, the solid lino shows the mean from our 10 Gaussian 
simulations and the dashed line is the corresponding analytic 
prediction of Kayo et al. (2001). 



that the PDF can be well described by a lognormal distri- 
bution. Here we show the mean PDF from our 10 Gaussian 
simulations as the solid line, and the errors correspond to 
the estimated variance. The dashed line shows the analytic 
relation given by Kayo et al. (2001), which is shown to be 
accurate at 5 ^ —0.5. It can be seen that our simulation 
matches well the analytic prediction. The main point of the 
plot is to show that our mass resolution is good enough to 
properly resolve the over-dense tail of the probability distri- 
bution function, for a typical smoothing scale used later on 
in our analysis. 

Our main motivation is to see whether the type of non- 
Gaussian models considered here, ie. fjvL models, can pro- 
duce more massive structures than does the standard Gaus- 
sian model. In Fig. [2] we present the evolution of the dark 
matter halo mass function as a function of redshift. The 
upper panel shows the results with WMAP-1 constrained 
fjvL model parameters. A distinctive feature is that a posi- 
tive fjvL, which gives a more strongly disturbed gravitational 
field, results in the formation of increased numbers of viri- 
alized structures. This effect is also much stronger at higher 
redshifts, as expected from the theoretical models of Matar- 
rese et al. (2000). In the lower panel of Fig. O we show the 
ratio of halo mass functions of non-Gaussian, dn/dM{NG), 
and Gaussian models, dn/dM{G), by the thick lines with 
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Figure 2. Upper panel: dark matter halo mass functions with 
/jVL parameters 134 (red dashed), (black solid), -58 (blue 
dotted). Lower panel: the ratio of halo mass functions of non- 
Gaussian and Gaussian models. Thick lines with error bars 
are for our suite of DM simulations and thin lines are the 
analytic model predictions of Matarrese et al. (2000), for the 
same set of /nl parameters (same line styles and colours as 
in the upper panel). 



error bars from our ten statistically equivalent realizations. 
Even though the simulation predicts the ratio of halo mass 
function to increase with halo mass and with increasing red- 
shift, there is still a large discrepancy with the analytical 
model expectations of Matarrese et al. (2001). Our simula- 
tions produce a much larger number density of dark mat- 
ter halos for non-Gaussian models than expected by Press- 
Schechter-based model predictions. The main difference be- 
tween the Press-Schechter halo mass functions and the N- 
body mass functions can be attributed to a different treat- 
ment of the dynamics of halo formation (cf. Komatsu et al. 
2003). It is valid to ask whether primordial non-Gaussianity 
also changes the dynamics of halo formation. This issue is 
unfortunately beyond our current study, as following the 
dynamics of halo formation requires much higher-resolution 
simulations than the ones we have run, together with more 
regularly spaced simulation outputs. We therefore delay this 
aspect of non-Gaussianity to future work. 

3.2 Protoclusters at z 2-4 

As previously shown, there are significant differences be- 
tween the Gaussian model and the /jv_t models at high red- 
shift, in particular for large positive /jvl values, which re- 
sults in more structure being formed. Miley et al. (2004) have 
reported a protocluster at z = 4.1 around a Ly-a emitter 
galaxy. They found 21 spectroscopically confirmed galaxies 
around the Ly-Q emitter galaxy, and the 1-dimensional ve- 
locity dispersion among these galaxies is about 325 s~^ km. 
Kurk et al. (2004) also observed a protocluster ai z — 2.16 
with a velocity dispersion of ~ 360 s~^ km over 3 to 5 Mpc, 
corresponding to co-moving sphere radii of 5 to 9 h~^ Mpc. 
Following Robinson & Baker (2000) and Mathis, Diego & 
Silk (2004) , we sample such structures in our simulations by 
throwing down a large number of spheres with radii of 7 h~^ 
Mpc and calculate their 1-D dark matter velocity dispersion. 

In Fig. |3] we show the cumulative probability distribu- 
tion function of the measured 1-D dark matter velocity dis- 
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Figure 3. Cumulative probability distribution function of the 
one-dimensional velocity dispersion measured within spheres 
of radii 7 h~^ Mpc. Upper (lower) panel show the z = 4: {z = 
2) results. Vertical dashed lines indicate the observed velocity 
dispersion by Miley et al. (2004) and Kurk et al. (2004). 



persion. In the top panel, we compare our findings with the 
measurement of Miley et al. (2004) at z ~ 4 . In the bot- 
tom panel, we plot the equivalent curves ai z — 2 and com- 
pare with the velocity dispersion measurement of Kurk et al. 
(2004). Contrary to Mathis et al. (2004), we find that our 
Gaussian simulation with as ~ 0.9 can explain the existence 
of structures with velocity dispersions larger than 320 s~^ 
km, such as the protocluster of Miley et al. (2004) . Although 
the simulation of Mathis et al. (2004) used the same cosmo- 
logical parameters as our simulation, they adopt too small 
a simulation volume (box of 100 h~^ Mpc aside), with the 
consequence that the cut-off in the input power spectrum on 
large scales depresses the velocity dispersion in their simu- 
lation. 

Fig. [3] shows that for all simulations with ag — 0.9, it 
is possible to observe a z ~ 4 protocluster with 1-D ve- 
locity dispersion of about 320 s~^ km. On the other hand 
in a low as universe (like one with as ~ 0.74, as favoured 
by WMAP-3), the chance of observing such protoclusters is 
very unlikely, unless there is a strongly disturbed primor- 
dial non-Gaussian field (such as one with /nl = 134). We 
note that these are not in any case very strong constraints, 
as we have ignored any velocity bias between the galaxies 
and the dark matter, which to some extent could solve the 
apparent discrepancy. Therefore we have not attempted to 
quantify the predicted number density of such objects at 
z ~ 4. Only detection of a large number of protoclusters 
of similar size or bigger at z ~ 4 could actually put real 
constraints on this simple model, together with some theo- 
retical constraints on the strength of velocity bias, as would 
be measurable in higher resolution N-body simulations. The 
lower panel of Fig. [3] shows that, at 2 ~ 2, it is rather easy 
to find, in all of the simulations considered, regions with 
velocity dispersion larger than 400 s~^ km. This highlights 
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Figure 4. Top panel: Probability distribution function of the 
mass (galaxy) density contrast within spheres of 12 Mpc 
radii: four sets of dark matter PDFs from Gaussian and non- 
Gaussian models (see labels in panel) and, in magenta, the 
PDF of 2dFGRS galaxies. Lower panel: the ratio of the 
PDFs shown in the top panel to the Gaussian PDF. This ex- 
plicitly shows that the PDF of 2dFGRS galaxies is dis- 
tinctively different to any of the DM models presented and 
highlights particularly well that a lower erg value increases the 
discrepancy between data and models even further. 



the fact that only very high redshift protoclusters are likely 
to set constraints on the cosmological parameters using the 
1-D velocity dispersion. 



3.3 Large scale structures in the 2dFGRS 

A clear limitation of the preceding results is the sparseness 
and quality of the data that we are comparing to. Hence 
we consider in this section constraints from local large-scale 
structure results. We estimate the galaxy probability distri- 
bution function (PDF) using a 2dFGRS LJJ^ volume-limited 
sample (Norberg et al. 2002; Baugh et al. 2004), compara- 
ble in volume to a simulation cube of ~ 200 h^^ Mpc on 
the side. The PDF is measured by Fast-Fourier- Transform 
(FFT) smoothing of the galaxy distribution with a spherical 
top-hat filter of radius R. For the smoothing scales consid- 
ered (4 to 16 h~^ Mpc), we recover the results of Croton 
et al. (2004), who used a count-in-cells (CiC) method to 
measure the galaxy PDF. The FFT-based method, once we 
accurately account for edge effects (particularly important 
when smoothing over large scales), is several orders of mag- 
nitude faster than standard CiC. 

In the top panel of Fig. |4l we show the 2dFGRS L^j 
galaxy PDF, smoothed over 12 h~^ Mpc, together with the 
dark matter PDF of the /ivL models specified in the caption. 
In the bottom panel, we show the ratio between the PDFs 
presented in the top panel and the PDF of the Gaussian 
simulation. Clearly the 2dFGRS galaxy PDF lies well above 
all model predictions, and in particular it is very discrepant 
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Figure 5. Top panel: Probability distribution function of the 
mass (galaxy) density contrast within spheres of 12 h~^ Mpc 
radii: in black, the PDF of a /jvl = 500 model; in blue, a 
biased ACDM PDF; in magenta, the PDF of 2dFGRS L* 
galaxies. Both models assume rrg = 0.90 cosmology. Lower 
panel: the ratio of the PDFs shown in the top panel to the 
Gaussian ACDM PDF. For a high erg value, both models seem 
able to reproduce the tail of the data PDF rather well. 



with models using a low ag value, as currently favoured by 
WMAP-3 data. 

Two possible approaches may help reconcile the models 
with the observations. Firstly, in Fig. |4l we compare the 
galaxy PDF with the mass PDF, ignoring any bias between 
them. We consider a simple linear bias model (b ~ 1.2) 
between galaxy and dark matter on the weakly non-linear 
scales of 12 h~^ Mpc in the Gaussian model. Secondly, we 
consider a more extreme positive non-Gaussian model with 
/jvL = 500, which predicts a significantly higher probability 
in the tail of the PDF. The constraints from WMAP on 
/jvL are only for very large scales, and the non-Gaussianity 
may be scale-dependent. In Fig. (5] we compare these two 
alternative models with the 2dFGRS L^^ PDF. 

In the following discussion, we assume that the data 
extracted from the 2dFGRS is a fair representation of the 
Universe, i.e. we require our model to easily reproduce the 
observations, something one potentially expects from a vol- 
ume close to 8 X 10** h"^ Mpc^. Clearly, if ag = 0.9, then 
both a standard ACDM with relatively strong linear bias 
on non-linear scales {b ~ 1.2) or a strong non-Gaussian 
model (/jvL = 500) without bias reproduce the 2dFGRS 
PDF rather well. In the case of a low as Universe, the 
bias needed for a pure Gaussian model would be uncom- 
fortably large (6 ~ 1.46), especially when accounting for the 
most recent determination of galaxy bias on weakly non- 
linear scales using the same 2dFGRS samples (Gaztanaga 
et al. 2005). Alternatively, a strong positive non-Gaussianity 
(./jvL ~ 500) coupled with a moderate galaxy bias could re- 
solve the discrepancy between data and models. In fact, a 
linear bias model on the scales of 12 h~^ Mpc may not be 
realistic. Gaztanaga et al. (2005) have recently measured the 
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galaxy bias on weakly non-linear scales (ie. ~ 10 to ~ 30 h^^ 
Mpc) using the 2dFGRS samples and found that the bias 
is non-linear in over-dense regions {5 S> 1). However the 
second order corrections are found to be negative, lowering 
the galaxy density, and hence favouring the non-gaussian 
interpretation. 

Two issues are worth mentioning in this analysis; we 
focus in particular on the LJJj volume-limited sample, be- 
cause of the presence of two very large coherent superstruc- 
tures in each 2dFGRS regionsQ, which strongly influences all 
higher-order statistics (Baugh et al. 2004; Croton et al. 2004; 
Gaztanaga et al. 2005; Nichol et al. 2006), and even to some 
extent the 2-point correlation function. The other reason is 
that Lj^ J galaxies are expected to only be very mildly biased 
tracers of the underlying dark matter (eg. Verde et al. 2002; 
Gaztanaga et al. 2005). Results from a larger volume-limited 
sample such as the ~ 2.5 L^^ sample are more difficult to 
interpret, mostly because the galaxy density is much lower, 
making the tail of the PDFs harder to measure. 



show how difficult it is to reproduce the observed high tail 
in the galaxy PDF, when smoothed on 12 h~^ Mpc scales. 
In a (T8 = 0.90 Universe, we argued that only a Gaussian 
ACDM model with strong linear bias between galaxies and 
mass (6 ~ 1.2) or a very strongly non-Gaussian model, with 
/jvL ~ 500, can reconcile the results mapped by the 2dF- 
GRS. The bias constraints obtained by Gaztanaga et al. 

(2005) using the same 2dFGRS samples favour the latter 
option. However the complication of galaxy bias on weakly 
non-linear scales in over-dense regions {5 ^ 1) makes it 
difficult to constrain the model parameters based on a sim- 
ple linear bias model. Nevertheless, in a low erg Universe as 
favoured by recent WMAP-3 results, the case for primordial 
non-Gaussianity is strengthened. Further tests, with real- 
istic galaxy formation recipes, are needed to properly ad- 
dress these issues. Our conclusions complement those in a 
recent theoretical paper by Hikage, Komatsu & Matsubara 

(2006) , namely that non-Gaussianity can be well constrained 
by high redshift galaxy surveys and that non-linear galaxy 
bias makes it difficult to detect primordial non-Gaussianity 
at low redshift. 



4 CONCLUSIONS 

In this paper, we have studied a simple non-Gaussian model 
in which the initial potential has a small perturbation rela- 
tive to the Gaussian random field and can be described by a 
simple non-linear coupling parameter (Komatsu et al. 2003) . 
We produce a series of N-body simulations to model struc- 
ture formation with both Gaussian and non-Gaussian initial 
conditions. We estimate the predicted halo mass function as 
a function of redshift, and make comparison with observa- 
tions for the PDF of the 1-D DM velocity dispersion as well 
as for the PDF of the DM overdensity. 

We find that, compared to the Gaussian model, a non- 
Gaussian model with deeper potential wells (ie. positive 
/jvl) produces larger structures, whereas models with neg- 
ative /nl values do the opposite. The differences are very 
pronounced at high redshift, but are less evident at z ~ 0. 
The presence of large coherent structures at high redshift is 
a better indicator of early structure formation than similar 
objects at low redshift. Under the assumption of no veloc- 
ity bias between galaxy and mass, our results show that it 
would be difficult to observe the z ^ 4.1 protocluster of 
Miley et al. (2004) in a erg = 0.74 Universe with Gaussian 
initial conditions, but becomes much easier with a strong 
non-Gaussian field (eg. /ml ~ 134). Protoclusters at high 
redshift can put strong constraints on non-Gaussian mod- 
els, once the efi'ects of bias and sample selections are prop- 
erly addressed. High redshift structure formation potentially 
provides a more powerful test of possible primordial non- 
Gaussianity than does the CMB, albeit on smaller scales. 

This of course assumes that we live in a typical region 
of the universe. However the large-scale coherent structures 
observed in the 2DFGRS (Baugh et al. 2005) that domi- 
nate the higher order statistics, and also seen in the SDSS, 
might indicate that a larger volume of the universe must be 
sampled to robustly test non-Gaussianity. 

A detailed comparison between the 2dFGRS L^^ galaxy 
PDF and the PDFs from a series of dark matter models 
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